
function [avg_diff_N,avg_diff_O,wage_sample_N1,wage_sample_N2,wage_sample_O1,wage_sample_O2,pr,pr_qrtr,pr_noncollege,pr_college,mm,wp_results,wp_noncollege_results,wp_college_results,wg_results,wg_noncollege,wg_college] = DO_moments_func(h_grid,simul,simul_size,param,dist_U,dist_N,dist_O,opt_ind_hp_N,opt_ind_hp_O,opt_p_U,opt_p_N,opt_p_O)
    
    %----- Specify Parameters -----%
    nu = param(1);
    lambda = param(6);          % Probability of search OTJ
    delta =param(11);           % Separation rate

    nphi=2;
     % Specify grids %
    hlb = 0.0;                 % Human capital grid lower bound 
    hub = 7.0;                 % Human capital grid upper bound 
    nh = 3000;                 % Number of possible human capital values
    h_grid = linspace(hlb,hub,nh); % Create human capital grid
    nt = 400;               % Number of quarters to follow each individual - 100 yrs max - most will exit the model before this
    %% =============================================================================================================== %%
    %========================================== Compute Moments of Interest ============================================%
    %================================================================================================================= %%
    
    % ----------- m1: JTJ Probability
    jtj_rate = 0.0;
    m_sum = 0.0;
    for ind_phi=1:nphi
        for ind_h = 1:nh
    
            ind_hp_N = opt_ind_hp_N(ind_phi,ind_h); 
            ind_hp_O = opt_ind_hp_O(ind_phi,ind_h);
                
            jtj_rate = jtj_rate + lambda*opt_p_N(ind_phi,ind_hp_N)*(1.0-delta)*dist_N(ind_phi,ind_h)*(1.0-nu);
            jtj_rate = jtj_rate + lambda*opt_p_O(ind_phi,ind_hp_O)*(1.0-delta)*dist_O(ind_phi,ind_h)*(1.0-nu);
    
            jtj_rate = jtj_rate + opt_p_U(ind_phi,ind_hp_N)*delta*dist_N(ind_phi,ind_h)*(1.0-nu);
            jtj_rate = jtj_rate + opt_p_U(ind_phi,ind_hp_O)*delta*dist_O(ind_phi,ind_h)*(1.0-nu);
            
            m_sum = m_sum + dist_N(ind_phi,ind_h)*(1.0-nu) + dist_O(ind_phi,ind_h)*(1.0-nu);
        end 
    end 
       
    if (m_sum>0.0) 
        jtj_rate = (jtj_rate/m_sum);
    else
        jtj_rate = 999.0;
    end
    m1 = jtj_rate;
    
    % ----------- m2: Unemployment Rate
    m2 = (sum(sum(dist_U(:,:)))/( sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:))) + sum(sum(dist_U(:,:))) ))*100.0;
    
    % ----------- m3: Total Percent Outsourced
    m3 = (( sum(sum(dist_O(:,:))) )/(sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:)))))*100.0;
    
    % ----------- m4: Ratio Outsourced With to Without Bachelor's Degree
    DO_phi1 = (( sum(dist_O(1,:)) )/( sum(dist_O(1,:)) + sum(dist_N(1,:)) ))*100.0;
    DO_phi2 = (( sum(dist_O(2,:)) )/( sum(dist_O(2,:)) + sum(dist_N(2,:)) ))*100.0;
    m4 = DO_phi1/DO_phi2;
    
    % ----------- m5: Avg. Quarterly Separation Rate
    delta_rate = 0.0;
    m_sum = 0.0;
    
    for ind_phi=1:nphi
        for ind_h = 1:nh
    
            ind_hp_N = opt_ind_hp_N(ind_phi,ind_h); 
            ind_hp_O = opt_ind_hp_O(ind_phi,ind_h);    
    
            delta_rate = delta_rate + (1.0-opt_p_U(ind_phi,ind_hp_N))*delta*dist_N(ind_phi,ind_h)*(1.0-nu);
            delta_rate = delta_rate + (1.0-opt_p_U(ind_phi,ind_hp_O))*delta*dist_O(ind_phi,ind_h)*(1.0-nu);
    
            m_sum = m_sum + dist_N(ind_phi,ind_h)*(1.0-nu) + dist_O(ind_phi,ind_h)*(1.0-nu);
        end 
    end 
    m5 = delta_rate/m_sum;
    
    % ----------- m6: Average Wage Loss After Unemployment Spell   
    wg_uspell = 0;
    m_sum = 0;
    for id = 1:simul_size
        stat_u = 0;
        for t = 1:nt-2 % Need to observe employed wage, unemployment, and then employed wage again
            
            if (simul(id,t).stat~=999)
                if (stat_u==1 && simul(id,t).stat~=0) % Was Unemployed but found Employment
                    w_after = simul(id,t).wage;
                    wg_uspell = wg_uspell + (log(w_after) - log(w_before))*100;
                    m_sum = m_sum +1;
                    stat_u = 0;
                end
    
                if (simul(id,t).stat~=0 && simul(id,t+1).stat==0) % Employed in t and unemployed in t+1
                    w_before = simul(id,t).wage;
                    stat_u = 1;
                end 
            end % End if has not exited the model (if stat~=999)
        end
    end
    wg_uspell = wg_uspell/m_sum;
    m6 = wg_uspell;
    
    % ----------- m7: Wage Dispersion p90/p50
    count_wages = 0;
    for id = 1:simul_size
        for t = 1:nt       
            if (simul(id,t).stat~=0 && simul(id,t).stat~=999)
                count_wages = count_wages + 1;
            end    
        end
    end

    % --- First save sample of wages
    wage_sim = zeros(1,count_wages);
    count_wages = 0;
    for id = 1:simul_size
        for t = 1:nt       
            if (simul(id,t).stat~=0 && simul(id,t).stat~=999)
                count_wages = count_wages + 1;
                wage_sim(count_wages) = simul(id,t).wage;
            end    
        end
    end

    p90_wage = prctile(wage_sim, 90);
    p50_wage = prctile(wage_sim, 50);
  
    m7 = p90_wage/p50_wage;
    
    
    %----------- m8: Average Annual Real Wage Growth
    % average wage growth over 1 year
    wg_sum2 = 0.0;
    m_sum = 0.0;
    for id = 1:simul_size
        for t = 1:(nt-4)
            if (simul(id,t).stat~=0 && simul(id,t+4).stat~=0 && simul(id,t).stat~=999 && simul(id,t+4).stat~=999) 
                w1 = simul(id,t).wage;
                w5 = simul(id,t+4).wage; 
                wg_sum2 = wg_sum2 + (log(w5) - log(w1))*100.0;
                m_sum = m_sum + 1.0;
            end 
        end 
    end 
    wg_sum2 = wg_sum2/m_sum;
    m8 = wg_sum2;
    
    % ----------- m9: College Wage Premium
    avg_wage_phi1 = 0.0;
    avg_wage_phi2 = 0.0;
    sum_phi1= 0.0;
    sum_phi2= 0.0;
    
    for id = 1:simul_size
        for t = 21:nt  % Age >= 25 yrs
            
            if (simul(id,t).stat~=0 && simul(id,t).stat~=999)   % If Not Unemployed or dead
                if (simul(id,t).phi_ind == 1) 
                    avg_wage_phi1 = avg_wage_phi1 + simul(id,t).wage;
                    sum_phi1 = sum_phi1 + 1.0;
                else
                    avg_wage_phi2 = avg_wage_phi2 + simul(id,t).wage;
                    sum_phi2 = sum_phi2 + 1.0;
                end 
            end 
        end 
    end 
    
    avg_wage_phi1 = avg_wage_phi1/sum_phi1;
    avg_wage_phi2 = avg_wage_phi2/sum_phi2;
    
    m9 = avg_wage_phi2/avg_wage_phi1;
    
    %---------- m10 - Unemployment Rate Among Age 20
    count_E =0.0;
    count_U =0.0;
    for id = 1:simul_size
        t=1; % age 20
        if (simul(id,t).stat~=999)
            if (simul(id,t).stat==0)  
                count_U = count_U + 1.0;
            elseif (simul(id,t).stat==1 || simul(id,t).stat==2) 
                count_E = count_E + 1.0;
            end    
        end 
    end
    
    m10 = (count_U/(count_E + count_U))*100.0;
    
    
    mm =[m1;m2;m3;m4;m5;m6;m7;m8;m9;m10]; % Model Moments
    %------------ Untargeted: Wage Penalty and Relative Wage Growth -----------------%
    
    count_total = 0;
    count_wp_noncollege =0;
    count_wp_college = 0;
    
    count_growth = 0;
    count_growth_noncollege=0;
    count_growth_college=0;
    
    t_save = zeros(id,1);
    for id = 1:simul_size
        % Randomly choose t to observe respondent
        r = rand;
        t = round(r*nt);
        t_save(id) = t;
        if (t>0)
            if (simul(id,t).stat~=999 && simul(id,t).stat~=0) % if employed
                count_total = count_total + 1;
                if (simul(id,t).phi_ind==1)
                    count_wp_noncollege = count_wp_noncollege + 1;
                else
                    count_wp_college = count_wp_college + 1;
                end
            end
        end
        if (t>0 && t+4<=nt)
            if (simul(id,t).stat~=999 && simul(id,t).stat~=0 && simul(id,t+4).stat~=999 && simul(id,t+4).stat~=0) % if employed now and next yr
                count_growth = count_growth + 1;
                if (simul(id,t).phi_ind==1)
                    count_growth_noncollege = count_growth_noncollege + 1;
                else
                    count_growth_college = count_growth_college + 1;
                end
            end
        end
    end % end loop over id
    
    X_wp = zeros(count_total,7);
    X_wp_noncollege = zeros(count_wp_noncollege,6);
    X_wp_college = zeros(count_wp_college,6);
    
    wage_vect = zeros(count_total,1);
    wage_vect_noncollege = zeros(count_wp_noncollege,1);
    wage_vect_college = zeros(count_wp_college,1);
    
    X_wg = zeros(count_growth,10);
    X_wg_noncollege = zeros(count_growth_noncollege,9);
    X_wg_college = zeros(count_growth_college,9);
    
    wg_vect = zeros(count_growth,1);
    wg_vect_noncollege = zeros(count_growth_noncollege,1);
    wg_vect_college = zeros(count_growth_college,1);
    
    count_total=0;
    count_wp_noncollege =0;
    count_wp_college = 0;
    count_growth=0;
    count_growth_noncollege=0;
    count_growth_college = 0;
    
    for id = 1:simul_size
        t = t_save(id);
    
        if (t>0)
            if (simul(id,t).stat~=999 && simul(id,t).stat~=0)
                count_total=count_total+1;
                
                wage_vect(count_total,1) = log(simul(id,t).wage);
                if  (simul(id,t).stat==2)
                    X_wp(count_total,1) = 1;
                end
                if  (simul(id,t).stat==1)
                    X_wp(count_total,1) = 0;
                end
                
                if simul(id,t).phi_ind==2
                    X_wp(count_total,2) = 1;                          % college
                else
                    X_wp(count_total,2) = 0;                          % non-college
                end       
                years = 20 + (t-1)/4;
                X_wp(count_total,3) = years-18;                       % exp_yr
                X_wp(count_total,4) = (years-18)^2;                   % exp_yr_sq
                X_wp(count_total,5) = h_grid(simul(id,t).h_ind);      % h
                X_wp(count_total,6) = log(h_grid(simul(id,t).h_ind)); % lnh
                X_wp(count_total,7) = 1; %constant
                
                % Save Results for Non-College Type
                if (simul(id,t).phi_ind==1)
                    count_wp_noncollege = count_wp_noncollege + 1;
                    wage_vect_noncollege(count_wp_noncollege,1) = log(simul(id,t).wage);
                    % Indicator for outsourced or not
                    if  (simul(id,t).stat==2)
                        X_wp_noncollege(count_wp_noncollege,1) = 1;
                    end
                    if  (simul(id,t).stat==1)
                        X_wp_noncollege(count_wp_noncollege,1) = 0;
                    end
                    % Controls
                    years = 20 + (t-1)/4;
                    X_wp_noncollege(count_wp_noncollege,2) = years-18;                       % exp_yr
                    X_wp_noncollege(count_wp_noncollege,3) = (years-18)^2;                   % exp_yr_sq
                    X_wp_noncollege(count_wp_noncollege,4) = h_grid(simul(id,t).h_ind);      % h
                    X_wp_noncollege(count_wp_noncollege,5) = log(h_grid(simul(id,t).h_ind)); % lnh
                    X_wp_noncollege(count_wp_noncollege,6) = 1;                              % constant
                end
                
                % Save Results for College Type
                if (simul(id,t).phi_ind==2)
                    count_wp_college = count_wp_college + 1;
                    wage_vect_college(count_wp_college,1) = log(simul(id,t).wage);
                    % Indicator for outsourced or not
                    if  (simul(id,t).stat==2)
                        X_wp_college(count_wp_college,1) = 1;
                    end
                    if  (simul(id,t).stat==1)
                        X_wp_college(count_wp_college,1) = 0;
                    end
                    % Controls
                    years = 20 + (t-1)/4;
                    X_wp_college(count_wp_college,2) = years-18;                           % exp_yr
                    X_wp_college(count_wp_college,3) = (years-18)^2;                       % exp_yr_sq
                    X_wp_college(count_wp_college,4) = h_grid(simul(id,t).h_ind);          % h
                    X_wp_college(count_wp_college,5) = log(h_grid(simul(id,t).h_ind));     % lnh
                    X_wp_college(count_wp_college,6) = 1;                                  % constant
                end
            end
        end
    
        if (t>0 && t+4<=nt)
            if (simul(id,t).stat~=999 && simul(id,t).stat~=0 && simul(id,t+4).stat~=999 && simul(id,t+4).stat~=0) % if employed
                count_growth = count_growth + 1;
                
                wg_vect(count_growth,1) = log(simul(id,t+4).wage) - log(simul(id,t).wage);
                % Indicator for outsourced or non-outsourced for both observations %
                if  (simul(id,t).stat==2 && simul(id,t+4).stat==2)
                    X_wg(count_growth,1) = 1;
                end
                if  (simul(id,t).stat==1 && simul(id,t+4).stat==1)
                    X_wg(count_growth,1) = 0;
                end
                % Indicator for move from outsourced to non-outsourced
                if  (simul(id,t).stat==2 && simul(id,t+4).stat==1) 
                    X_wg(count_growth,2) = 1; 
                else
                    X_wg(count_growth,2) = 0;
                end
                % Indicator for college or non-college
                if simul(id,t).phi_ind==2
                    X_wg(count_growth,3) = 1;                            % college
                else
                    X_wg(count_growth,3) = 0;                            % non-college
                end
                % Other controls
                years = 20 + (t-1)/4;
                X_wg(count_growth,4) = years-18;                         % exp_yr
                X_wg(count_growth,5) = (years-18)^2;                     % exp_yr_sq
                X_wg(count_growth,6) = h_grid(simul(id,t+4).h_ind);      % h
                X_wg(count_growth,7) = log(h_grid(simul(id,t+4).h_ind)); % lnh
                X_wg(count_growth,8) = simul(id,t).wage;
                X_wg(count_growth,9) = (simul(id,t).wage)^2;
                X_wg(count_growth,10) = 1;                               %constant
                
                
                if (simul(id,t).phi_ind==1) %----- Save data for non-college group only
                    count_growth_noncollege = count_growth_noncollege + 1;
                    wg_vect_noncollege(count_growth_noncollege,1) = log(simul(id,t+4).wage) - log(simul(id,t).wage);
                    % Indicator for outsourced or non-outsourced for both observations %
                    if  (simul(id,t).stat==2 && simul(id,t+4).stat==2)
                        X_wg_noncollege(count_growth_noncollege,1) = 1;
                    end
                    if  (simul(id,t).stat==1 && simul(id,t+4).stat==1)
                        X_wg_noncollege(count_growth_noncollege,1) = 0;
                    end
                    % Indicator for move from outsourced to non-outsourced
                    if  (simul(id,t).stat==2 && simul(id,t+4).stat==1) 
                        X_wg_noncollege(count_growth_noncollege,2) = 1; 
                    else
                        X_wg_noncollege(count_growth_noncollege,2) = 0;
                    end
                    % Other controls
                    years = 20 + (t-1)/4;
                    X_wg_noncollege(count_growth_noncollege,3) = years-18;                         % exp_yr
                    X_wg_noncollege(count_growth_noncollege,4) = (years-18)^2;                     % exp_yr_sq
                    X_wg_noncollege(count_growth_noncollege,5) = h_grid(simul(id,t+4).h_ind);      % h
                    X_wg_noncollege(count_growth_noncollege,6) = log(h_grid(simul(id,t+4).h_ind)); % lnh
                    X_wg_noncollege(count_growth_noncollege,7) = simul(id,t).wage;
                    X_wg_noncollege(count_growth_noncollege,8) = (simul(id,t).wage)^2;
                    X_wg_noncollege(count_growth_noncollege,9) = 1;                               %constant
                    
                else                        %----- Save data for college group only
                    count_growth_college = count_growth_college + 1;
                    wg_vect_college(count_growth_college,1) = log(simul(id,t+4).wage) - log(simul(id,t).wage);
                    % Indicator for outsourced or non-outsourced for both observations %
                    if  (simul(id,t).stat==2 && simul(id,t+4).stat==2)
                        X_wg_college(count_growth_college,1) = 1;
                    end
                    if  (simul(id,t).stat==1 && simul(id,t+4).stat==1)
                        X_wg_college(count_growth_college,1) = 0;
                    end
                    % Indicator for move from outsourced to non-outsourced
                    if  (simul(id,t).stat==2 && simul(id,t+4).stat==1) 
                        X_wg_college(count_growth_college,2) = 1; 
                    else
                        X_wg_college(count_growth_college,2) = 0;
                    end
                    % Other controls
                    years = 20 + (t-1)/4;
                    X_wg_college(count_growth_college,3) = years-18;                         % exp_yr
                    X_wg_college(count_growth_college,4) = (years-18)^2;                     % exp_yr_sq
                    X_wg_college(count_growth_college,5) = h_grid(simul(id,t+4).h_ind);      % h
                    X_wg_college(count_growth_college,6) = log(h_grid(simul(id,t+4).h_ind)); % lnh
                    X_wg_college(count_growth_college,7) = simul(id,t).wage;
                    X_wg_college(count_growth_college,8) = (simul(id,t).wage)^2;
                    X_wg_college(count_growth_college,9) = 1;    
                end
                
            end
        end
    end
    % Wage penalty regression -- No Sample Restrictions
    wp_results = regress(wage_vect,X_wp);
    
    % Wage penalty regression -- Only Non-College
    wp_noncollege_results = regress(wage_vect_noncollege,X_wp_noncollege);
    
    % Wage penalty regression -- Only College
    wp_college_results = regress(wage_vect_college,X_wp_college);
    
    % Wage growth regression -- No Sample Restrictions
    wg_results = regress(wg_vect,X_wg);
    
    % Wage growth regression -- Only Non-College
    wg_noncollege = regress(wg_vect_noncollege,X_wg_noncollege);
    
    % Wage growth regression -- Only College
    wg_college = regress(wg_vect_college,X_wg_college);


    %=================================================================%
    %===== Use simul to compute annual transition probabilities ======%
    %=================================================================%
    sum_Nu = 0;
    sum_Ou = 0;
    sum_NO = 0;
    sum_ON = 0;
    sum_uN = 0;
    sum_uO = 0;

    sum_u = 0; % count those starting in unemployment
    sum_N = 0; % count those starting non-outsourced
    sum_O = 0; % count those starting outsourced

    for id = 1:simul_size
        for t = 1:(nt-4)    % Need to look 4 quarters ahead

            if ((simul(id,t).stat==0) && (simul(id,t+4).stat~=999))     % Was unemployed in first year observed and had not exited the model by the next year
                sum_u = sum_u + 1;
                if (simul(id,t+4).stat==1)            % Was non-outsourced one year later
                    sum_uN = sum_uN + 1;
                elseif (simul(id,t+4).stat==2)        % Was outsourced one year later
                    sum_uO = sum_uO + 1;
                else
                end

            elseif ((simul(id,t).stat==1) && (simul(id,t+4).stat~=999)) % Was non-outsourced in the first year observed and had not exited the model by the next year
                sum_N = sum_N + 1;
                if (simul(id,t+4).stat==0)            % Was unemployed one year later
                    sum_Nu = sum_Nu + 1;
                elseif (simul(id,t+4).stat==2)        % Was outsourced one year later
                    sum_NO = sum_NO + 1;
                else
                end

            elseif ((simul(id,t).stat==2) && (simul(id,t+4).stat~=999)) %Was outsourced in the first year observed and had not exited the model by the next year
                sum_O = sum_O + 1;
                if (simul(id,t+4).stat==0)            % Was unemployed one year later
                    sum_Ou = sum_Ou + 1;
                elseif (simul(id,t+4).stat==1)        % Was non-outsourced one year later
                    sum_ON = sum_ON + 1;
                else
                end
           
            else % stat==999 (dead)
            end

        end
    end

    pr(1) = sum_Nu/sum_N;
    pr(2) = sum_Ou/sum_O;

    pr(3) = sum_NO/sum_N;   
    pr(4) = sum_ON/sum_O;

    pr(5) = sum_uN/sum_u;
    pr(6) = sum_uO/sum_u;

    %----------- Transition Probabilities - Only Non-college
    sum_Nu = 0;
    sum_Ou = 0;
    sum_NO = 0;
    sum_ON = 0;
    sum_uN = 0;
    sum_uO = 0;

    sum_u = 0; % count those starting in unemployment
    sum_N = 0; % count those starting non-outsourced
    sum_O = 0; % count those starting outsourced

    for id = 1:simul_size
        for t = 1:(nt-4)    % Need to look 4 quarters ahead
            if (simul(id,t).phi_ind==1) % Only Non-College
                if ((simul(id,t).stat==0) && (simul(id,t+4).stat~=999))     % Was unemployed in first year observed and had not exited the model by the next year
                    sum_u = sum_u + 1;
                    if (simul(id,t+4).stat==1)            % Was non-outsourced one year later
                        sum_uN = sum_uN + 1;
                    elseif (simul(id,t+4).stat==2)        % Was outsourced one year later
                        sum_uO = sum_uO + 1;
                    else
                    end

                elseif ((simul(id,t).stat==1) && (simul(id,t+4).stat~=999)) % Was non-outsourced in the first year observed and had not exited the model by the next year
                    sum_N = sum_N + 1;
                    if (simul(id,t+4).stat==0)            % Was unemployed one year later
                        sum_Nu = sum_Nu + 1;
                    elseif (simul(id,t+4).stat==2)        % Was outsourced one year later
                        sum_NO = sum_NO + 1;
                    else
                    end

                elseif ((simul(id,t).stat==2) && (simul(id,t+4).stat~=999)) %Was outsourced in the first year observed and had not exited the model by the next year
                    sum_O = sum_O + 1;
                    if (simul(id,t+4).stat==0)            % Was unemployed one year later
                        sum_Ou = sum_Ou + 1;
                    elseif (simul(id,t+4).stat==1)        % Was non-outsourced one year later
                        sum_ON = sum_ON + 1;
                    else
                    end
           
                else % stat==999 (dead)
                end
            end
        end
    end

    pr_noncollege(1) = sum_Nu/sum_N;
    pr_noncollege(2) = sum_Ou/sum_O;

    pr_noncollege(3) = sum_NO/sum_N;   
    pr_noncollege(4) = sum_ON/sum_O;

    pr_noncollege(5) = sum_uN/sum_u;
    pr_noncollege(6) = sum_uO/sum_u;
    
    
    %----------- Transition Probabilities - Only College
    sum_Nu = 0;
    sum_Ou = 0;
    sum_NO = 0;
    sum_ON = 0;
    sum_uN = 0;
    sum_uO = 0;

    sum_u = 0; % count those starting in unemployment
    sum_N = 0; % count those starting non-outsourced
    sum_O = 0; % count those starting outsourced

    for id = 1:simul_size
        for t = 1:(nt-4)    % Need to look 4 quarters ahead
            if (simul(id,t).phi_ind==2) % Only College
                if ((simul(id,t).stat==0) && (simul(id,t+4).stat~=999))     % Was unemployed in first year observed and had not exited the model by the next year
                    sum_u = sum_u + 1;
                    if (simul(id,t+4).stat==1)            % Was non-outsourced one year later
                        sum_uN = sum_uN + 1;
                    elseif (simul(id,t+4).stat==2)        % Was outsourced one year later
                        sum_uO = sum_uO + 1;
                    else
                    end

                elseif ((simul(id,t).stat==1) && (simul(id,t+4).stat~=999)) % Was non-outsourced in the first year observed and had not exited the model by the next year
                    sum_N = sum_N + 1;
                    if (simul(id,t+4).stat==0)            % Was unemployed one year later
                        sum_Nu = sum_Nu + 1;
                    elseif (simul(id,t+4).stat==2)        % Was outsourced one year later
                        sum_NO = sum_NO + 1;
                    else
                    end

                elseif ((simul(id,t).stat==2) && (simul(id,t+4).stat~=999)) %Was outsourced in the first year observed and had not exited the model by the next year
                    sum_O = sum_O + 1;
                    if (simul(id,t+4).stat==0)            % Was unemployed one year later
                        sum_Ou = sum_Ou + 1;
                    elseif (simul(id,t+4).stat==1)        % Was non-outsourced one year later
                        sum_ON = sum_ON + 1;
                    else
                    end
           
                else % stat==999 (dead)
                end
            end
        end
    end

    pr_college(1) = sum_Nu/sum_N;
    pr_college(2) = sum_Ou/sum_O;

    pr_college(3) = sum_NO/sum_N;   
    pr_college(4) = sum_ON/sum_O;

    pr_college(5) = sum_uN/sum_u;
    pr_college(6) = sum_uO/sum_u;

    %=================================================================%
    %==== Use simul to compute QUARTERLY transition probabilities ====%
    %=================================================================%
    sum_Nu = 0;
    sum_Ou = 0;
    sum_NO = 0;
    sum_ON = 0;
    sum_uN = 0;
    sum_uO = 0;

    sum_u = 0; % count those starting in unemployment
    sum_N = 0; % count those starting non-outsourced
    sum_O = 0; % count those starting outsourced

    % Also compute average overall hc for next comparison for low vs high hc individuals %
    sum_hc = 0;
    count_hc = 0;
    for id = 1:simul_size
        for t = 1:(nt-1)    % Need to look 1 quarter ahead
            if (simul(id,t).stat~=999)
                sum_hc = sum_hc + h_grid(simul(id,t).h_ind);
                count_hc = count_hc + 1;
            end

            if ((simul(id,t).stat==0) && (simul(id,t+1).stat~=999))     % Was unemployed in first quarter observed and had not exited the model by the next quarter
                   
                sum_u = sum_u + 1;
                if (simul(id,t+1).stat==1)            % Was non-outsourced one quarter later
                    sum_uN = sum_uN + 1;
                elseif (simul(id,t+1).stat==2)        % Was outsourced one quarter later
                    sum_uO = sum_uO + 1;
                else
                end

            elseif ((simul(id,t).stat==1) && (simul(id,t+1).stat~=999)) % Was non-outsourced in the first quarter observed and had not exited the model by the next quarter
                sum_N = sum_N + 1;
                if (simul(id,t+1).stat==0)            % Was unemployed one quarter later
                    sum_Nu = sum_Nu + 1;
                elseif (simul(id,t+1).stat==2)        % Was outsourced one quarter later
                    sum_NO = sum_NO + 1;
                else
                end

            elseif ((simul(id,t).stat==2) && (simul(id,t+1).stat~=999)) %Was outsourced in the first quarter observed and had not exited the model by the next quarter
                sum_O = sum_O + 1;
                if (simul(id,t+1).stat==0)            % Was unemployed one quarter later
                    sum_Ou = sum_Ou + 1;
                elseif (simul(id,t+1).stat==1)        % Was non-outsourced one quarter later
                    sum_ON = sum_ON + 1;
                else
                end
           
            else % stat==999 (dead)
            end

        end
    end

    pr_qrtr(1) = sum_Nu/sum_N;
    pr_qrtr(2) = sum_Ou/sum_O;

    pr_qrtr(3) = sum_NO/sum_N;   
    pr_qrtr(4) = sum_ON/sum_O;

    pr_qrtr(5) = sum_uN/sum_u;
    pr_qrtr(6) = sum_uO/sum_u;
    avg_sim_hc = sum_hc/count_hc;

    %======================================================================%
    %==== Use simul to compute low vs high HC transition probabilities ====%
    %======================================================================%
    sum_uE_low = 0;
    sum_uE_high = 0;

    sum_u_low = 0;    % count those starting in unemployment
    sum_u_high = 0;

    for id = 1:simul_size
        for t = 1:(nt-1)    % Need to look 1 quarter ahead

            if ((simul(id,t).stat==0) && (simul(id,t+1).stat~=999))     % Was unemployed in first quarter observed and had not exited the model by the next quarter
                
                if ( h_grid(simul(id,t).h_ind) < avg_sim_hc)
                    sum_u_low = sum_u_low + 1;

                    if (simul(id,t+1).stat==1 || simul(id,t+1).stat==2 )     % Entered employment
                        sum_uE_low = sum_uE_low + 1;
                    end
                else
                    sum_u_high = sum_u_high + 1;

                    if (simul(id,t+1).stat==1 || simul(id,t+1).stat==2 )     % Entered employment
                        sum_uE_high = sum_uE_high + 1;
                    end
                end
           
            end
        end
    end

    tr_UE_low = sum_uE_low/sum_u_low;
    tr_UE_high = sum_uE_high/sum_u_high;



    %=================================================================%
    %======= For Figure 2: Save Normalized Wage Disrtibutions ========%
    %=================================================================%
    % Save wage samples for 4 different groups
    count_N1 = 0;
    count_O1 = 0;
    count_N2 = 0;
    count_O2 = 0;
    count_wages = 0;
    
    for id = 1:simul_size
        for t = 1:nt
            if (simul(id,t).stat~=999 && simul(id,t).stat~=0) % If employed
                count_wages = count_wages + 1;
    
                if (simul(id,t).phi_ind==1)    % Non-College
                    if (simul(id,t).stat==1)     % If Non-Outsourced
                        count_N1 = count_N1 + 1;
                    elseif (simul(id,t).stat==2) % If Outsourced
                        count_O1 = count_O1 + 1;
                    else
                    end                
                else                            % College
                    if (simul(id,t).stat==1)     % If Non-Outsourced
                        count_N2 = count_N2 +1;
                    elseif (simul(id,t).stat==2) % If Outsourced
                        count_O2 = count_O2 + 1;
                    else
                    end
                end
            end
        end % end loop over nt
    end % end loop over ids in simul_size
    
    
    wage_sample_total = zeros(1,count_wages);
    wage_sample_N1 = zeros(1,count_N1);
    wage_sample_O1 = zeros(1,count_O1);
    wage_sample_N2 = zeros(1,count_N2);
    wage_sample_O2 = zeros(1,count_O2);
    
    count_N1 = 0;
    count_O1 = 0;
    count_N2 = 0;
    count_O2 = 0;
    count_wages = 0;
    for id = 1:simul_size
        for t = 1:nt
            if (simul(id,t).stat~=999 && simul(id,t).stat~=0) % If employed
                count_wages = count_wages + 1;
                wage_sample_total(count_wages) = simul(id,t).wage;
                if (simul(id,t).phi_ind==1)      % Non-College
                    if (simul(id,t).stat==1)     % If Non-Outsourced
                        count_N1 = count_N1 + 1;
                        wage_sample_N1(count_N1) = simul(id,t).wage;
                    elseif (simul(id,t).stat==2) % If Outsourced
                        count_O1 = count_O1 + 1;
                        wage_sample_O1(count_O1) = simul(id,t).wage;
                    else
                    end                
                else                             % College
                    if (simul(id,t).stat==1)     % If Non-Outsourced
                        count_N2 = count_N2 +1;
                        wage_sample_N2(count_N2) = simul(id,t).wage;
                    elseif (simul(id,t).stat==2) % If Outsourced
                        count_O2 = count_O2 + 1;
                        wage_sample_O2(count_O2) = simul(id,t).wage;
                    else
                    end
                end
            end
        end % end loop over nt
    end % end loop over ids in simul_size
    
    avg_wage_total = mean(wage_sample_total); 
    
    wage_sample_N1(:) = wage_sample_N1(:)/avg_wage_total;
    wage_sample_N2(:) = wage_sample_N2(:)/avg_wage_total;
    wage_sample_O1(:) = wage_sample_O1(:)/avg_wage_total;
    wage_sample_O2(:) = wage_sample_O2(:)/avg_wage_total;


    %=================================================================%
    %========= For Figure 5: Save HC Investment Differences ==========%
    %=================================================================%
    opt_hp_N = h_grid(opt_ind_hp_N);
    opt_hp_O = h_grid(opt_ind_hp_O);
    h_diff_N(1,:) = opt_hp_N(1,:) - h_grid(:)';
    h_diff_N(2,:) = opt_hp_N(2,:) - h_grid(:)';

    h_diff_O(1,:) = opt_hp_O(1,:) - h_grid(:)';
    h_diff_O(2,:) = opt_hp_O(2,:) - h_grid(:)';

    % N vs O - College and Non-College
    avg_diff_N = zeros(2,12);
    avg_diff_O = zeros(2,12);
    for ind_phi = 1:nphi
        avg_diff_N(ind_phi,1) = mean(h_diff_N(ind_phi,1:108));        %0-0.25
        avg_diff_N(ind_phi,2) = mean(h_diff_N(ind_phi,109:215));      %0.25-0.5
        avg_diff_N(ind_phi,3) = mean(h_diff_N(ind_phi,216:322));      %0.5-0.75
        avg_diff_N(ind_phi,4) = mean(h_diff_N(ind_phi,323:429));      %0.75-1
    
        avg_diff_N(ind_phi,5) = mean(h_diff_N(ind_phi,430:536));      %1-1.25
        avg_diff_N(ind_phi,6) = mean(h_diff_N(ind_phi,537:643));      %1.25-1.5
        avg_diff_N(ind_phi,7) = mean(h_diff_N(ind_phi,644:750));      %1.5-1.75
        avg_diff_N(ind_phi,8) = mean(h_diff_N(ind_phi,751:857));      %1.75-2
    
        avg_diff_N(ind_phi,9) = mean(h_diff_N(ind_phi,858:964));      %2-2.25
        avg_diff_N(ind_phi,10) = mean(h_diff_N(ind_phi,965:1072));    %2.25-2.5
        avg_diff_N(ind_phi,11) = mean(h_diff_N(ind_phi,1073:1179));   %2.5-2.75
        avg_diff_N(ind_phi,12) = mean(h_diff_N(ind_phi,1180:1286));   %2.75-3
    
        avg_diff_O(ind_phi,1) = mean(h_diff_O(ind_phi,1:108));        %0-0.25
        avg_diff_O(ind_phi,2) = mean(h_diff_O(ind_phi,109:215));      %0.25-0.5
        avg_diff_O(ind_phi,3) = mean(h_diff_O(ind_phi,216:322));      %0.5-0.75
        avg_diff_O(ind_phi,4) = mean(h_diff_O(ind_phi,323:429));      %0.75-1
    
        avg_diff_O(ind_phi,5) = mean(h_diff_O(ind_phi,430:536));      %1-1.25
        avg_diff_O(ind_phi,6) = mean(h_diff_O(ind_phi,537:643));      %1.25-1.5
        avg_diff_O(ind_phi,7) = mean(h_diff_O(ind_phi,644:750));      %1.5-1.75
        avg_diff_O(ind_phi,8) = mean(h_diff_O(ind_phi,751:857));      %1.75-2
    
        avg_diff_O(ind_phi,9) = mean(h_diff_O(ind_phi,858:964));      %2-2.25
        avg_diff_O(ind_phi,10) = mean(h_diff_O(ind_phi,965:1072));    %2.25-2.5
        avg_diff_O(ind_phi,11) = mean(h_diff_O(ind_phi,1073:1179));   %2.5-2.75
        avg_diff_O(ind_phi,12) = mean(h_diff_O(ind_phi,1180:1286));   %2.75-3
    
    end


end


